//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// cl_fit_mv - Ajusta curva de luz de binarias a dados observados utilizando
//			   Cadeias de Markov. 
//
// Arquivo com a definicao dos cabecados das funcoes basicas.
// Definicao das funcoes em:
//
// - cl_fit_fit_funcs.cpp
// - cl_fit_input_funcs.cpp
// - cl_fit_data_funcs.cpp
// 
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

#include <iostream>
#include <cmath>
#include <fstream>
#include <cmath>
#include <vector>
#include <sstream>
#include <string>
#include <cstring>
#include <iomanip>
#include <gsl/gsl_siman.h>
#include <gsl/gsl_multimin.h>

//#include "cl_fit_input_funcs.cpp"
//#include "cl_fit_data_funcs.cpp"

//#include "light_curve.h"
//#include "light_curve.cpp"

//#include "cl_fit_mc_v0_1.hpp"
//#include "lib_zabot/dataFileReader.h"
#include "dataFileReader_v0.1.1/dataFileReader.h"

//----------------------------------------------------------------------------------------------------------------//
//
//----------------------------------------------------------------------------------------------------------------//
//
//----------------------------------------------------------------------------------------------------------------//
// ifndef serve para evitar que cabecalho seja definido mais de uma vez depois da chama da funcao...

#ifndef FIT_MC_HEADER

#define FIT_MC_HEADER

using namespace std;

void help_page(int level);

struct Data {
	double fase;
	double flux;
	double ferr;
};

/*
struct DataKeeper {
	vector<vector<double> > par;
	//  vector<double> vfac;
	vector<vector<Data> > cur_luz;
	//  vector<double> fase;
	//  vector<double> flux;
	//  vector<double> ferr;
	double f2max;
	ofstream* saida;
};
*/

//----------------------------------------------------------------------------------------------------------------//
//
//----------------------------------------------------------------------------------------------------------------//
//
//----------------------------------------------------------------------------------------------------------------//

class DataReader {
	
	int _nData; // numero de curvas de luz de entrada
	int* _npData; // numero de pontos por curva de luz
	
	string _ilist;
	string* _filename;
	Data** _input_; // Guarda todas as curvas de luz
	
public:
    DataReader(char* argv);
	
	int nData(void){return _nData;}; // Devolve numero de curvas de lux
	int nData(const int ncl){return _npData[ncl];}; // Devolce numero de pontos em cada curva de luz
	
	double inline DataReaderFase(const int i, const int j){return _input_[i][j].fase;};
	double inline DataReaderFlux(const int i, const int j){return _input_[i][j].flux;};
	double inline DataReaderFErr(const int i, const int j){return _input_[i][j].ferr;};
	
	void info(void);
};


//----------------------------------------------------------------------------------------------------------------//
//
//----------------------------------------------------------------------------------------------------------------//
//
//----------------------------------------------------------------------------------------------------------------//


class PReader {
	
	// Classe para ler parametros de entrada
	
	//Algumas definicoes do arquivo de entrada
	// Input File Definition - IFD:  
	//     Number of fixed Parameters lines      -  IFD_NP
	//     Number of columns Parameters          - *IFD_P (array[IFD_NP])
	//     Spot Line(onde com. def. de spot - 1) -  IFD_SL
	//     Number of Spot Parameters             -  IFD_SP
	//     posicao da lei de obs. de Borda       -  IFD_B
	
	int __IFDNP;        //= 2;
	int *__IFDP; //= {5, 7};
	int __IFDSL;        //= 3;
	int __IFDSP;        //= 5;
	int *__IFDB_2;         //= 2;
	int *__IFDB_1;
	double _f0;
	double _f1;
	double _df;
	vector<vector<double> > _par;
	vector<vector<string> > _pardef;
	vector<double> phases;
public:
    
	PReader(char* argv);
	
	int inline ifdnp(void){return __IFDNP;};
	int inline ifdp(const int i){return __IFDP[i];};
	int inline ifdb_1(const int i){return __IFDB_1[i];};
	int inline ifdb_2(const int i){return __IFDB_2[i];};
	int inline ifdsp(void){return __IFDSP;};
	
	string inline pardef(const int i, const int j){return _pardef[i][j];};
	
    int inline slaces(void){return int(_par[0][0]);};
    int inline definicoes_obscurecimento(const int law);    
    
	vector<vector<double> > inline org_par(int component);
	
    int nspot(int component = 0){return int(_par[2+component][3]);};
	//int nspot_2(void){return int(_par[3][3]);};
    double Iwd(void){return _par[3][0];};
    double Irr(int component = 0){return _par[2+component][7];};
    double Albedo(int component = 0){return _par[2+component][6];};
    double nfases(void){return phases.size();};
    double ff(int nf){return phases[nf];};
    double f0(void){return _f0;};
    double f1(void){return _f1;};
    double df(void){return _df;};
    void f0(double val){_f0=val;};
    void f1(double val){_f1=val;};
    void df(double val){_df=val;};
	
	double inline getpar(const int i,const int j){return _par[i][j];};
	void inline setpar(const int i,const int j,const double nval){_par[i][j] = nval;};
	
    double fot_erro(void)
    {
		if(_par[1][4] > 0.0) return (_par[2][0]+_par[2][7])/_par[1][4];
		else return 0.0;
    };
	
    void inline show_par(ostream& out)
    {
		for( int i = 0; i < _par.size(); i++){
			out<<"# ";
			for(int ii = 0; ii < _par[i].size(); ii++){
				if(i < 6) out << _pardef[i][ii] << _par[i][ii] << " | ";
				else out << _pardef[6][ii] << _par[i][ii] << " | ";
				//if(i == 4 || i == 5) out << _pardef[i][ii] << _par[i][ii] << " | ";
				
			}
			out<<endl;
		}
    };
};


//----------------------------------------------------------------------------------------------------------------//
//
//----------------------------------------------------------------------------------------------------------------//
//
//----------------------------------------------------------------------------------------------------------------//


class SimAnSolve : public DataReader, public PReader {

	//----------------------------------------------------------------------------------------------------//
	// Definicoes para gerador de numeros randomicos
	
	const gsl_rng_type * _Trng;
	gsl_rng * _rand;

	//----------------------------------------------------------------------------------------------------//
	// Determina quais parametros serao ajustados. Em principio, todos os parametros podem ser ajustados. Mas isso eh
	// altamente nao recomendavel. Na definicao da funcao irei criar um default para os principais parametros e isso
	// podera ser auterado apos a definicao da funcao utilizando FitPar(&)
	
	double** _FitPar; // quais parametros ajustar
	int** _FitPar_index; // indices para os parametros a ajustar
	int _NFitPar; // numero de parametros a ajustar
	
	double** _PrivParams; // Armazena parametros privados de cada curva de luz
	int _NPrivPar; // Numero de parametros privados
	int* _PrivPar_index; // posicao de cada parametro no vetor de parametro geral
		
	
	//----------------------------------------------------------------------------------------------------//
	// Parametros para o ajuste
	
	int _nPerturb; // numero maximo de perturbacoes por iteracao 
	int _nSuccess; // numero de sucessos por iteracao (diminui temperatura apos atingir esse numero de sucessos)
	int _nIter; // numero de iteracoes antes de finalizar programa. As interacoes sao contadas a cada troca de temperatura.
	int _nTrys; // numero de perturbacoes antes de avaliar eficiencia da cadeia. (~ 1/10 de _nPerturb)
	int _iter; // conta o numero de iteracoes
	
	double _T0; // Temperatura inicial da simulacao
	double _Tf; // Temperatura final da simulacao. Este parametro eh redundante com _nIter.
				// O programa atualiza um de acordo com o outro
	double _dT; // Fator multiplicativo (<1) para atualizar temperatura
	double _defautSSize; // default and uniform step size
	double* _StepSize; // quantos 1/100 eh o tamanho do passo em relacao a cada parametro (1 = 100 por cento)
	double* _MaxStepSize; // Maior passo que pode ser ajustado. 2x tamanho de passo default e eh diminuido a cada iteracao
	double* _MinMaxStepSize; // Paro de ajustar tamanho maximo do passo quando chego neste tamanho 1/10 do tamanho default

	
	//----------------------------------------------------------------------------------------------------//
	// Guarda informacoes de iteracoes
	
	string _OutputRootFilename;	// raiz do arquivo para guardar as solucoes
					// programa ira guardar solucao em "'root'_'itr'.dat" onde itr eh a iteracao
					// por default a raiz eh = 'simansolve'
	double** _Solutions;	// Armazera infos sobre cada iteracao em temperatura. Resultado eh armazenado 
							// em arquivo ao final de uma iteracao e vetor eh zerado.
	
public:
	SimAnSolve(char* arg1, char* arg2); //:PReader(arg1),DataReader(arg2);
	
	//
	//----------------------------------------------------------------------------------------------------//
	//----------------------------------------------------------------------------------------------------//
	// Definindo recuperacao e atualizaco de variaveis internas
	//
	
	int inline nFitPar(void){return _NFitPar;};
	
	int inline nTrys(void){return _nTrys;};
	void inline nTrys(const int nval){_nTrys = nval;};
	
	int inline NPerturb(void){return _nPerturb;}; // devolve numero de perturbacoes 
	void inline NPerturb(const int nval)
	{
		if( nval > NSuccess() ) _nPerturb = nval;
		else cout << "# - NPerturb[" << nval << "] must be larger than Nsuccess["<< NSuccess() <<"]. Cannot change value ...\n";
			
	};	// muda numero de perturbacoes
	
	int inline NSuccess(void){return _nSuccess;}; // devolve numero de sucessos
	void inline NSuccess(const int nval)
	{
		if( nval < NPerturb() ) _nSuccess = nval;
		else cout << "# - NPerturb[" << NPerturb() << "] must be larger than Nsuccess["<< nval <<"]. Cannot change value ...\n";

	};	// muda numero de sucessos

	int inline iter(void){return _iter;}; // devolve numero de iteracoes
										  // Vou utilizar o operador ++ para gerar uma iteracao
	
	int inline nIter(void){return _nIter;}; // devolve numero TOTAL de iteracoes
	void inline nIter(const int nval)
	{
		_nIter = nval;
		_Tf = _T0 * pow(dT(), double(nval));
		
	};	// muda numero TOTAL de iteracoes e atualiza Tf
	
	
	double inline T0(void){return _T0;}; // devolve Temperatura inicial 
	void inline T0(const double nval)
	{
		if (_iter == 0) {
			_T0 = nval;
			_Tf = nval * pow(dT(), double(nIter())); // atualiza Tf
		}
		else cout << "# [WARNING] - Iterations started. 'Public' function cannot change T0 ...\n";
		

		
	};	// muda Temperatura inicial. Funciona apenas antes de iniciar iteracoes!
	
	double inline dT(void){return _dT;}; // devolve fator de atualizacao da temperatura
	void inline dT(const double nval)
	{
		if (0.0 < nval && nval < 1.0) {
			_dT = nval;
			_Tf = T0() * pow(nval, double(nIter())); // atualiza Tf
		}
		else {
			cout << "# [WARNING] - 0 < nval[" << nval << "] < 1 - Cannot change value ...\n";
		}

	}; // muda fator de atualizacao da temperatura
	
	void inline Set_Priv_Par(int i, int j, double val){_PrivParams[i][j] = val;}; // muda valor do parametro privado
	double inline Get_Priv_Par(int i, int j){return _PrivParams[i][j];}; // acessa valor do parametro privado

	
	//
	//----------------------------------------------------------------------------------------------------//
	//----------------------------------------------------------------------------------------------------//
	// Definindo funcoes de vizualizacao e controle 
	//
	
	void inline show_FitPar(ostream& out);
	void inline show_FitPar2(ostream& out);
	void inline FitPar(const int npar, int* i_pos,int* j_pos);

	//
	//----------------------------------------------------------------------------------------------------//
	//----------------------------------------------------------------------------------------------------//
	// Definindo funcoes de minimizacao
	
	int inline Perturba(void); // Gera perturbacao nos parametros de entrada
	double inline Evaluate(void); // calcula o valor da funcao custo (dados - modelo) -> Chi^2
	int inline SelectPar(void);
	double inline dealer(const int index,double val,double step);
	
	int operator++ ();	// Pre incremento. Vou utilizar esse operador para gerar uma iteracao.
						// Por iteracao eu considero o processo a mesma temperatura, o que significa um laco de nSuccess(<200) ou 
						// nPerturb (<1000). Durante esse processo ele testa se os passos estao muito grande ou nao 
						// (a cada 100 perturbacoes).
						// Devolve _iter, o numero de iteracoes realizadas.
	
};

//----------------------------------------------------------------------------------------------------------------//
//
//----------------------------------------------------------------------------------------------------------------//
//
//----------------------------------------------------------------------------------------------------------------//

int definicoes_obscurecimento(const int law);
void show_par(ostream& out, vector<vector<double> > &par);
vector<vector<double> > org_par(vector<vector<double> > &in_par); //organiza os parametros para 
//ficarem de acordo requisitado pela classe Light_Curve.

void show_par2(ostream& out, vector<vector<double> > &par);

double qui2(void *);
void do_step(const gsl_rng *, void *, double);
double how_far(void *, void *);
void print_cfg(void *);
void copy(void *source, void *copy);
void* create_cfg(void *xp);
void destroy(void *xp);
void siman_solve(const gsl_rng *, void *, double (*)(void *), void (*)(const gsl_rng *, void *, double), double (*)(void *, void *), void (*)(void *), gsl_siman_params_t, void (*)(void *,void *),void* (*)(void *), void (*)(void *));
template <class T> double oculta(vector<double>& pos, vector<double>& fase, T& STAR);

template <class T>
double oculta(vector<double>& pos, vector<double>& fase, T& STAR)
{
	return STAR.oculta(pos,fase);
}

//----------------------------------------------------------------------------------------------------------------//
//
//----------------------------------------------------------------------------------------------------------------//
//
//----------------------------------------------------------------------------------------------------------------//


class GSL_Simplex : public DataReader, public PReader {

//----------------------------------------------------------------------------------------------------------------//
// Determina quais parametros serao ajustados. Em principio, todos os parametros podem ser ajustados. Mas isso eh
// altamente nao recomendavel. Na definicao da funcao irei criar um default para os principais parametros e isso
// podera ser auterado apos a definicao da funcao utilizando FitPar(&)
	
	double** _FitPar; // quais parametros ajustar
	int** _FitPar_index; // indices para os parametros a ajustar
	int _NFitPar; // numero de parametros a ajustar
	
	double** _PrivParams; // Armazena parametros privados de cada curva de luz
	int _NPrivPar; // Numero de parametros privados
	int** _PrivPar_index; // posicao de cada parametro no vetor de parametro geral
	
//----------------------------------------------------------------------------------------------------------------//
// Definicoes para chamada de funcoes da biblioteca GSL para minimizacao pelo procedimento de simplex.
	
	const gsl_multimin_fminimizer_type *T;
	gsl_multimin_fminimizer *s;
	gsl_vector *ss, *x;
	gsl_multimin_function minex_func;

	//----------------------------------------------------------------------------------------------------//
	// Parametros para o ajuste
		
	double _defautSSize; // default and uniform step size
	double* _StepSize; // quantos 1/100 eh o tamanho do passo em relacao a cada parametro (1 = 100 por cento)
	size_t _iter,_nIter;
	int _status;

	string _OutputRootFilename;	// raiz do arquivo para guardar as solucoes
								// programa ira guardar solucao em "'root'_'itr'.dat" onde itr eh a iteracao
								// por default a raiz eh = 'simansolve'
	
public:
	GSL_Simplex(char* arg1, char* arg2); //:PReader(arg1),DataReader(arg2);
	
	//
	//----------------------------------------------------------------------------------------------------//
	//----------------------------------------------------------------------------------------------------//
	// Definindo funcoes de vizualizacao e controle 
	//
	
	void inline show_FitPar(ostream& out);
	void inline show_FitPar2(ostream& out);
	void inline FitPar(const int npar, int* i_pos,int* j_pos, const int npriv, int* i_priv,int* j_priv);
	
	//
	//----------------------------------------------------------------------------------------------------//
	//----------------------------------------------------------------------------------------------------//
	// Definindo funcoes de minimizacao
	
	double inline Evaluate(const gsl_vector *par); // calcula o valor da funcao custo (dados - modelo) -> Chi^2
	double inline dealer(const int index,double val);
	void inline init(void); // Inicializa minimizador
	int operator++ ();	// Pre incremento. Vou utilizar esse operador para gerar uma iteracao.
	int inline iter(void){return _iter;}; // devolve numero de iteracoes
										  // Vou utilizar o operador ++ para gerar uma iteracao
	int inline status(void){return _status;};
    double inline simplex_size(void){return gsl_multimin_fminimizer_size(s);};	
	double inline score(void){return s->fval;};
	
	void inline store_res(void);
	void inline store_cl(void);
};

#endif